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DESIGNING ILLUMINATION LENSES AND MIRRORS BY THE 
NUMERICAL SOLUTION OF MONGE AMPERE EQUATIONS 
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Abstract. We consider the inverse refractor and the inverse reflector problem. The task is to 
design a free-form lens or a free-form mirror that, when illuminated by a point light source, produces 
a given illumination pattern on a target. Both problems can be modeled by strongly nonlinear 
second-order partial differential equations of Monge-Ampere type. In [Math. Models Methods Appl. 
Sci. 25 (2015), pp. 803-837, DOI: 10.1142/S0218202515500190 the authors have proposed a B-spline 
collocation method which has been applied to the inverse reflector problem. Now this approach is 
extended to the inverse refractor problem. We explain in depth the collocation method and how to 
handle boundary conditions and constraints. The paper concludes with numerical results of refracting 
and reflecting optical surfaces and their verification via ray tracing. 
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1. Introduction. Both problems, the inverse refractor and the inverse reflector 
problem, from illumination optics can be formulated in the following framework: Let 
a point-shaped light source and a target area be given, e.g. a wall. Then we would 
like to construct an apparatus that projects a prescribed illumination pattern, e.g. an 
image or a logo, onto the target. Since we aim for maximizing the efficiency, we would 
like to construct our optical device in such a way that, neglecting losses, it redirects 
all light emitted by the light source to the target. We focus our attention to the design 
of such an optical system in the simple case that it either consists of a single free-form 
lens or of a single free-form mirror, see Figure O for an illustration of the former 
case. Our goal is now to compute the shape of the optically active surfaces, modeled 
as free-form surfaces, such that the desired light intensity distribution is generated 
on the target. Since these problems from illumination optics from the mathematical 
point of view conceptually fall into the class of inverse problems, they are also called 
inverse reflector problem and inverse refractor problem, respectively. In particular, 
since the size of the optical system is comparable to that of the projected image, we 
address the case of the near field problems. 

There is a variety of technical applications of such optical systems, e.g. spotlights 
with prescribed illumination patterns used in street lamps or car headlamps, see 
e.g. [T1I51I5S]. 

The authors present in a solution method for the inverse reflector problem 
via numerically solving a strongly nonlinear second-order partial differential equation 
(PDE) of Monge-Ampere type. Due to the high potential of this approach we now 
extend this method to the case of illumination lenses. 
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Figure 1.1: Setting of the refractor problem. The index of refraction of the lens ma¬ 
terial is ni, while the surrounding has the refractive index n 2 - 


This paper is organized as follows: Since the reflector problem has been discussed 
in detail in [5] we mainly focus on the refractor problem. We start with the state of 
the art for its solution in Section Then, we formulate the problem via a partial 
differential equation of Monge-Ampere type which we discuss in Section for the 
construction of a refractor. For completeness we also give the Monge-Ampere formu¬ 
lation for the reflector problem in Section]^ Next, the numerical method is explained 
in Section]^ Since this type of optical design problem raises many difficulties in the 
solution process we discuss in Section how these can be resolved. Finally, in Sec¬ 
tion we look at numerical results for the inverse reflector and refractor problems 
and end this paper in Section with our conclusions. 


2. State of the art. In this section we discuss the methods available for the 
solution of the inverse design problems in nonimaging optics, see the monographies by 
Chaves [5] and by Winston, Mifiano and Benitez [SH] for an introduction to nonimaging 
optics and the paper by Patow and Pueyo |30j for a survey article on inverse surface 
design from the graphics community’s point of view. For a detailed survey of solution 
techniques for the inverse reflector problem, we refer the reader to [SI Section 2]. 

Focusing on the inverse refractor problem, in the paper by Wester an Bauerle m 
there is a list of approaches, a discussion on practical problems, e.g. extended sources 
and Fresnel losses, and examples with LED lighting, e.g. a lens for automotive fog 
light and a lens producing a logo. 

In the rest of this section, we first give a short overview of other solution techniques 
in Section 12.11 and then focus on methods based on PDEs in Section 12.21 which is 


also our problem formulation of choice. Finally, we discuss some advanced topics in 
Section 113] and draw our conclusions in Section HaI 


2.1. Approaches for the solution of inverse problems in nonimaging 
optics not based on a PDE. We distinguish three different groups of techniques 
for the solution of inverse problems in nonimaging optics, which are not based on a 
PDE: there are methods resorting from optimization techniques, others built from 
Cartesian ovals and a third group of methods which are geometrical constructions. 

Optimization approaches. There are methods for the design of optical surfaces, 
which are based on optimization techniques, see e.g. [33 HOI- Starting from an initial 
guess, the outline of the iterative optimization process for the determination of the 
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optical surfaces is as follows: First, the current approximation of the optical surfaces 
is validated by ray tracing. In a second step, using an objective function, which is 
often closely related to the Euclidean norm, the resulting irradiance distribution is 
compared to the desired one and a correction of the optical surfaces is determined. 
The process ends, when a suitable quality criterion is fulfilled, otherwise these two 
steps are repeated. 

The advantage of this method is that it is very flexible. However, optimization 
procedures are very costly because of the repeated application of the ray tracing and 
it is unclear if the iterative methods converge at all. 

Cartesian ovals methods. Cartesian ovals are curves of fourth order in the plane. 
They can be associated with two foci such that light emitted at one focus is collected 
at the other focus. Here the Cartesian oval coincides with the interface of two optical 
media with different refractive indices. Cartesian ovals can be extended to surfaces in 
3d with the same property. By combining clippings of several of these surfaces in an 
iterative procedure a new segmented surface can be constructed that approximates 
the solution. This strategy has first been developed by Kochengin and Oliker [20l [21] 
for the construction of solutions for the inverse reflector problem. Later this has been 
extended to the inverse refractor problem using Cartesian ovals, see also [H 
Section 2] for some theoretical background, and for a collimated light beam instead 
of a point light source |5S| using hyperboloids. 

Although this technique has the advantage to permit the construction of contin¬ 
uous but non-differentiable surfaces [25| . the number of clippings K required grows 
linearly in the number of pixels in the image. For example, using ellipsoids of rev¬ 
olution for the construction of a mirror with accuracy 7 > 0 , the complexity of the 
method scales like 0{^ log ^), see [ 22 ], such that it quickly becomes infeasible for 
higher resolutions. 

Geometric construction methods. Reflective and refractive free-form surfaces can 
also be designed by geometric approaches. Probably the most famous of these tech¬ 
niques is the simultaneous multiple surfaces (SMS) method extending the ideas of 
Cartesian-oval methods, see e.g. [39l Chapter 8 ] and [3l |24| and the references therein. 
The main idea of the SMS method is the simultaneous construction of two optical sur¬ 
faces, e.g. both surfaces of a lens, which permits to couple two prescribed incoming 
wave fronts, e.g. coming from two point light sources, with two prescribed outgoing 
wave fronts. While in its 2D version, the method is used to design rotationally sym¬ 
metric optical surfaces, in a 3D variant it is also capable to construct free-form optical 
surfaces. However, the authors could not find any hint on the computational costs 
in the literature but conjecture that this scheme is expensive especially for complex 
target illumination patterns. 

2.2. Solution techniques via PDE approaches. In several publications for 
the inverse refractor problem a PDE is derived, whose solution models the desired 
optical free-form surfaces, see e.g. [Ml na miiiziiiii Ea ESI HD ■ In these approaches 
usually the low wavelength limit is assumed to hold, i.e. the problems are formulated 
using the geometrical optics approximation. 

Some examples for the inverse refractor problem with a more complex target 
illumination pattern are shown in [Ml SIl EH ES] • However, in all four articles the 
descriptions and discussions of the numerical methods are incomplete. To the best of 
the authors’ knowledge the solution method is not fully documented in the literature. 

While we consider the case of a point light source, an interesting and closely 
related problem is shaping the irradiance distribution of a collimated light beam, see 
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e.g. [MIES for the theory including some results on existence and uniqueness of 
solutions. 

We refer the reader to the monography by Gutierrez |14j for a general overview 
of Monge-Ampere-type equations. 

Since we are looking for an optical surface which redirects light coming from a 
source onto a target, one can model this problem in terms of optimal transportation. 

Optimal transport. There are also methods which are based on a problem of 
optimal transport which leads to Monge-Ampere-type equations, see e.g. UlliEH]. 
First the ray mapping, i.e. the mapping of the incoming light rays onto the points at 
the target, is computed via an optimal transport approach. At this point the optical 
surface is still unknown but in a next step it is constructed from the knowledge of the 
target coordinates for each incoming light ray. In 1998 Parkyn [53] already described 
a very similar procedure. 

2.3. Advanced topics. In the current formulation of the problem only one 
single idealized point light source has been used. An extension to multiple point light 
sources is discussed by Lin [23] where the optical refractors are determined from those 
calculated for single point light sources by a weighted least-squares approach. More 
techniques for the case of extended light sources can be found in the papers by Bortz 
and Shatz [3| and Wester et al. [55] . 

In particular for the refractor problem, some energy is lost for the illumination of 
the target because of internal reflections in the lens material. A theoretical discussion 
of these Fresnel losses can be found in the publications by Gutierrez m Section 5.13] 
and Gutierrez and Mawi |17j . In [TJIS] the losses are minimized by free-form shaping 
of both refractive surfaces of the lens. 


2.4. Conclusion. Our approach is motivated by the fact that even for the special 
case of a single point light source and the computation of just one surface of the lens we 
could not find any fully detailed method in the literature which can produce complex 
illumination patterns on the target area. From the authors’ point of view, the most 
promising approach is the one by solving a PDF of Monge-Ampere type. 


3. The inverse refractor problem. This section is devoted to the formulation 
of the Monge-Ampere equation that models the near field refractor problem as given 
in the paper by Gutierrez and Huang m- Since the full theory is a bit involved, we 
restrict ourselves to a summary of the most important aspects and refer the reader to 
m Appendix A] and the paper by Karakhanyan and Wang [19] for the details. Our 
notation also follows these sources. 

We now proceed as follows: At first, we fix the geometric setting and the implicit 
definition of the refracting and the target surfaces in Section[3T] Then we apply Snell’s 


law of refraction in Section 3.2 and follow the path of the light ray in Section 3.3 
Finally, in Section [3^ we obtain the desired equation of Monge-Ampere type. 


3.1. The Geometric Setting. Since a lens has two surfaces we need to design 
both of them. For simplicity we choose a spheric inner surface, i.e. the surface which 
faces the light source is a subset of a sphere with center at the position of the light 
source. Thus there is no refraction of the incoming light at this interface, the inner 
surface is optically inactive. 

It remains to compute the shape of the outer surface facing the target area. 
To that end let us define the quotient k = ^ of the refractive indices of the lens 
material ni and the environment n 2 . We assume that the light source illuminates a 
non-empty subset U of the northern hemisphere of the unit sphere C The 
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third component of a n incoming light ray with direction x = (cci,a; 2 ,G U is 
then given as X3 = Thus we define x' := (xi,X2)^ and parametrize 

our outer lens surface by the distance function p := p(x'), i.e. the surface is given as 
r := {p(x')x : X e U}. 

The target E is defined as a subset of a hypersurface implicitly given by the zero 
level set of a continuously differentiable function ip via 

(3.1) E C {z e . ^( 2 ;) = 0}. 

Note that for the numerical solution procedure in the Newton-type method we require 
that Ip is twice continuously differentiable. While in general much more complicated 
situations are supported [TB] , for simplicity we restrict ourselves to the case where the 
target E is on a shifted x-y-plane such that ip{z) := Z3 — 7 for a shift 7 > 0. 

To model the luminous intensity of the source we define the density function 
f : U ^ K'*', where K’*' := {x S M : x > 0}. The corresponding density function for 
the desired illumination pattern on the target E is denoted by 5 : E —>■ K’*'. Since we 
want to redirect all incoming light onto the target the density functions need to fulfill 
the energy conservation condition 


(3.2) 


fdS= / gdS. 


Note that for simplicity we neglect the loss of reflected light intensity. For a more 
complicated derivation of a Monge-Ampere-type equation for the refractor problem 
taking losses into account see m- 

3 . 2 . Snell’s law of refraction. According to Snell’s law of refraction in vec¬ 
torial notation (see e.g. [HI Chapter 4.4] or [HI Chapter 12]), the direction of the 
light ray after refraction at the point px is y = ^(x — <i>(x • v)v) G S^, where 
<i>(s) := s — + s^ — 1 and i/ is the outer unit normal on T defined as a func¬ 

tion on U. 

As detailed in [HI (2.15)], for the outer normal unit vector at x G (7 we find 


(3.3) 


-Vp(x') 4- x(p(x') -L Vp(x') ■ x') 
^p^(x') + lVp(x')12 - (Vp(x') • x')2 ’ 


where V/ denotes the gradient of a function / and Vp(x') := (Vp(x'),0) G K^. To 
ease notation, we define the utility function G which represents the denominator in 
(ra, i.e. 


(3.4) 


G(x',u,p) := -L 1 p 12 - (p • xOA 


3.3. Following the light ray. Next, we consider the line which contains the 
light ray after refraction, defined by the point px and the direction vector y. We now 
turn to finding the point z = {zi,Z 2 , where the refracted light ray hits the target 
E. In order to determine the third component Z 3 of z, we first define the utility point 
w = {wi,W 2 , 0 )^ as the intersection point of this line with the plane {x S : X 3 = 0 } 
which is given as 


w = p(x')x -L doy 
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for a do G R- For a proof of the existence of w see [Ml Appendix A.2]. Using (3.3), 
we confirm that 


(3.5) 


= A(x',p(x'), Vp(x'))V(p2), 


where the utility function F is given by 


(3.6) 


A(x',u,p) 


1 $(“/G(x'.ii,p)) 

2 -G'(x', M, p) + (m + P • x')$(“/G(x',ii.p)) ' 


After refraction at the point px the light ray hits the target E at point z = (zi, Z 2 ^ z^Y' 
given as 


z(x') = p(x')x + diy = p(x')x + t(w - p(x')x). 


From the third component we know that t = • 

We introduce the short notation x 0 y := xy”^. Let us define F := F{x', p) 

and denote its partial derivatives by Dx'A, DpF and DpF, respectively. 

A lengthy computation using standard calculus and some tensor identities of 
Sherman-Morrison type yields 

Dw' = 2pFMB^p + B 


where M := I + p®T)pF and B := 2FVp(8)Vp + V(p^)(8)Dx'A + DuFV(p^)(8)Vp. 
Note that 


M-^ =1- 


Vp 0 DpF 
F + Vp • DpF’ 


In a bit more involved computation along the same lines we compute 
Dz' = 2tpFM.{l — j3{W — px') • Vi/')(D^p + A) 


with /? := (V^ • (w — px)) := (V'pi: '*Ap 2 ) = -^(x^ ^p)j where 

A'.={tB + (1 — t)C) and 
2tpb 

C := D(px') H —(w' - px') (g) V(pa; 3 ). 

PX3 

3.4. Monge—Ampere equation. The energy conservation 
holds if we replace U with any arbitrary subset U C U and E with E := T{U) C E, 
where T : {7 —?► E, x i—?> z(x'). By coordinate transformation this yields the identity 
det(Dz) = f/{gY^ - |x'p). Finally, we can derive the Monge-Ampere equation for 
the refractor problem 


3.2) clearly also 


(3.7) 


det(Dp + A) 


/(x) 

p(z(x'))i7’ 


for x' S D 


where ft := {{xi,X 2 Y ^ • {xi,X 2 ,X 3 Y G U} H = id(x',p, Vp) is computed 

by 


H:={1- |x'n|VV'|(2t)2p3(-^)F(F + Vp • DpF), 


see [Ml Appendix A]. 
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Existence and uniqueness of solutions. In general, for boundary value problems 
with Monge-Ampere equations proving well-posedness, i.e. existence and uniqueness 
of the solution and continuous dependency on the parameters, is a hard problem, e.g. 
see [ini Section 1.4] for an example of a discretized Monge-Ampere equation obtained 
by finite differences on a grid of 4 x 4 cells which has 16 different solutions. 

Some theoretical results for the existence of a solution for the refractor problem 
under some appropriate conditions can be found in |16j in Theorem 5.8 for k < 1 and 
Theorem 6.9 for k > 1. Additionally there are results on the uniqueness of the solution 
if just finitely many single points on the target are illuminated, see [T6l Theorem 5.7] 
for K < 1 and m Theorem 6 . 8 ] for k > 1. 

For proving existence and uniqueness of a solution one typically requires the 
equation of Monge-Ampere type to be elliptic. A necessary condition is that the right- 
hand side of ( |3.7| ) is positive. For this reason we demand that /3 < 0 or, equivalently, 
VV’ • (w — px) < 0. If this term is positive we can simply replace "0 by —if. 

4. The inverse reflector problem. The inverse reflector problem can be mod¬ 
eled as a Monge-Ampere-type equation very similarly to the case of the inverse re¬ 
fractor problem in Section]^ see m- 

Using the same definitions and notation as in Section]^ and introducing the sub¬ 
stitution u := we first define 

p’ 


t := 1 — u —, 


M := I + 


X3 

X (X) X 


w := — Vu, and 
d 


d := |Vm1^ — {u — Vm • x)^ 
b := - (Vu-x)2, 

1 / 1 , 

z := —X -I- t Uw- X I . 

u \ u 


We assume that t > 0, i.e., ^ > zs, and VV'- (w — 7X) > 0 . Then the Monge-Ampere 
equation for the inverse reflector problem reads 


det 





(mw —x)-VV^ /(x) 

45 5(z) ’ 


see m and [S] for the details. 


5. Numerical solution of partial differential equations of Monge Am¬ 
pere type. The numerical solution of strongly nonlinear second-order PDFs, includ¬ 
ing those of Monge-Ampere type, is a highly active topic in current mathematical 
research. There are many different approaches available on the market, see the re¬ 
view paper by Feng, Glowinski and Neilan m and also [S] for an overview. 

However, most methods are not well-suited for all equations of Monge-Ampere 
type such that it remains unclear if a particular method can be successfully applied 
to our problems. In the authors propose to use a spline collocation method which 
turns out to provide an efficient solution strategy for Monge-Ampere equations arising 
in the inverse reflector problem. 


In Section |5.1| we explain the idea of a collocation method, which reduces the 
problem to finding an approximation of the solution within a finite dimensional space. 


Then we discuss the choice of appropriate basis functions in Section 5.2 
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5.1. Collocation method. As discretization tool for the Monge-Ampere equa¬ 

tions arising in the reflector and refractor problem, we propose a collocation method, 
see e.g. Bellomo et al. [5] for examples of collocation methods applied to nonlinear 
problems. Let the PDE F{x, u, Vu, D^m) = 0 in fl and constraints G(x, u, Vu) = 0 on 
dfl be given. In this setting we approximate rt in a finite-dimensional trial subspace 
of G^(r2), i.e. for some finite set I and basis functions C G^(n) we choose the 

ansatz ii = Next, we only require that the PDE holds true on a collocation 

set D C D which contains only finitely many points. So our approximation ii of the 
solution of our PDE satisfies 

F(r,M(T), Vw,D^'u) = 0 , for r e D, 

G(t, u(r), Vu) = 0, for r S i9D. 

This discrete nonlinear system of equations is solved by a quasi-Newton method, 
which uses trust-region techniques for ensuring global convergence of the method, see 
Chapter 4.2.1 in [5] and the references cited therein for the details and the proofs. 

5.2. Splines and collocation points. We choose to apply a space of spline 
functions as ansatz space because of their advantageous properties, see e.g. [TOll^ldK] 
for details on the theory of splines. 

For a given interval [a, 5] we fix an equidistant knot sequence T = with 

n-fold knots at the interval end points a = ti for 1 < f < n and b = ti for -|- 1 < 
i < N + n. Moreover, we require that the knot sequence is strictly increasing inside 
the interval (a, b), i.e. U < for n < i < N. 

Then, an appropriate basis for our spline space is given by the B-spline functions 
Ni^n of order n which can be defined via the recursion formula 

= X[U,U+^]{t), = (A’j.n-l * 

where X[ti,ti+i] is the characteristic function of the interval C K and the 

convolution of two functions is defined as (/ * g)(x) := J^f{s)g{x — s) ds. Since we 
require that the ansatz functions are twice differentiable, we choose cubic splines, i.e. 
n = 4. 

In two dimensions the ansatz functions on a rectangular domain are obtained via a 
tensor ansatz and then are used as the Bi in the previous subsection. The collocation 
points are chosen to coincide with the sequence of equidistant knots. Since this leads 
to an underdetermined system of equations we use a not-a-knot condition at the very 
but last knot at each interval end, i.e. we require that the spline function is three times 
continuously differentiable at this knot. In other words, the restriction of the spline 
to the union of the two subintervals closest to each interval end is a cubic polynomial 
and the knot could be removed without changing the spline function. This is a much 
simpler approach than the one used in the previous work 0 Section 4.2.3] but provides 
approximately the same accuracy. 

6 . Numerical solution of equations of Monge Ampere type for optical 
applications. Next, we consider the particular difficulties that we have to overcome 
to efficiently solve the equations of Monge-Ampere type that arise in the reflector and 
refractor problems. 

6.1. Boundary conditions. The boundary conditions for both, the inverse re¬ 
flector and refractor problems, are realized via a Picard-type iteration as similarly 
proposed by Froese m Section 3.4]. 
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We assume that the light rays hitting the boundary of the optical surface also 
hit the boundary of the target, i.e. z{dfl) = i9E, see Section 4.5] for the details. 
This assumption is related to the edge ray principle, see e.g. [321 Appendix Bj. Note 
that z also depends on the solution p and its derivative Vp. In order to have a 
boundary condition which is easier to handle, we do not fix the target coordinate 
on the boundary but only its normal component. Since we do not know the normal 
component of the mapping z for the exact solution we proceed as follows: For solving 
the nonlinear system of equations from our collocation method we use a Newton-type 
method producing iterations p^, k = 1, 2,..., Umax, starting with an initial guess p°. 
We denote the corresponding mappings by z^ := z(x', p^, Vp^). In the kth iteration 
we require that the outer normal of the mapping of the current iteration and of 
the orthogonal projection of the mapping z^~^ of the last iteration onto the boundary 
coincide, i.e. 


z — arg min z — z 
zeas 


fe-i I 


i/(x') = 0 


for x' € dfl, 


see m Section 4.5] (cf. also |T31 Section 3.3]). The left-hand side is then used as 
function G in Section [5Tl 

Since the last iteration is involved in the boundary condition the function G 
changes in each iteration so that we solve different problems in successive steps. In 
order to ensure the existence of a solution of the subproblems we follow the approach 
by Froese [H Section 3.4] and add a parameter c in front of the right-hand side of 
the Monge-Ampere equation (3.7), i.e. we replace / by cf where c is an additional 
unknown in our equation. An additional constraint to compensate this new degree of 
freedom is discussed in Section lOl 


6.2. Ellipticity constraint. For proofs of results for existence and uniqueness 
of a solution we require the equation of Monge-Ampere type to be elliptic. In order to 
ensure ellipticity we manipulate the determinant in the same way as explained in |S1 
Section 4.4] (cf. also [T31 Section 4.3]): Let W = [Wij]i<ij <2 G be a matrix. 
For a penalty parameter A > 0 we define the modified determinant 


detJW := maxjO, Wiq} maxjO, W 2 , 2 } — Wi ,2 

- A [(min{0,>Vi,i})^ + (min{0, >V 2 , 2 })^] 


which we use instead of the determinant in the left-hand side of the Monge-Ampere 
equation (3.7). For an elliptic solution of this equation the left-hand side is exactly 
the same for the determinant and the modified determinant, see [SI Lemma 4.2]. 
Furthermore, each non-elliptic solution of (3.7) is not a solution of this equation, 
when the determinant is replaced by the modified determinant. 


6.3. Choice of the “size” of the refractor. Up to now the refractor is at 
most uniquely determined up to its size. Therefore we define our initial guess Mq 
of the problem appropriately and search for a solution u of same size requiring that 
f^uds = uq ds holds true, see also [SJ Section 4.5]. 

Note that this condition is taken account of by the additional unknown c intro¬ 
duced in Section l6T] 


6.4. Total internal reflection. In case that k < 1 it is possible that a ray 
of light exceeds the critical angle and total internal reflection occurs, such that this 
light ray does not reach the target. Of course we know that this is not true for the 
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solution, because we require that all light rays hit the target. However during the 
iteration process of our nonlinear solver this phenomenon can appear. If this is the 
case the argument of the square root in the definition of $ in Section |3.2| is negative 
at this position. 

To overcome this instability we replace ‘I*(s) by its stabilized counterpart $(s) := 
s — A/max{0, — 1}. Then the situation of total internal reflection is treated 

like the case when the light ray hits the surface exactly at the critical angle. The 
refracted light ray most likely also misses the target and therefore this intermediate 
step cannot satisfy the Monge-Ampere equation (3.7) such that further iterations are 
performed. 

If total internal reflection doesn’t occur, which is the case we intend to have for 
our solution, we have $( 5 ) = $(s) and therefore obtain an equivalent problem. 


6.5. Nested iteration. The convergence of Newton-type methods sensitively 
depends on the choice of an initial guess that is close enough to the solution. We 
apply a nested iteration strategy in order to largely increase the stability of the solver 
but also in order to accelerate the solution procedure. We start with a coarse grid for 
the spline surface and a blurred version of the image for the illumination pattern. 

The blurring process is necessary because a coarse grid cannot produce a very 
detailed image on the target area. For this reason we convolve the image, which is 
given as a raster graphic in our case, with a discrete version of the standard mollifier 
function i^(x) := exp(“Y(i-|xp)) if |x| < 1 and zero otherwise, namely with ipn{i,j) := 
for n G N and indices i,j G Z for the pixel coordinates. 

If our grid has N x N nodes we alternately increase the resolution N of the grid 
and decrease the strength n of blurring, i.e. we solve the problem for different pairs 
of {N,n), see also [51 Sections 4.3 and 5.2]. 


6 .6. Initial guess. For the refractor problem we simply use the surface of a 
sphere with center at the position of the light source and a prescribed radius as initial 
guess. 

For the reflector problem we start with a reflective surface producing a homoge¬ 
neous illumination pattern on the target. We obtain this reflector by first using the 
method of supporting ellipsoids uniisi] and our collocation technique afterwards, see 
also 0 Section 5.2.3]. 


6.7. Minimal gray value. The density function g corresponds to the target 
illumination on E and is given by 8 bit digital grayscale images (integer gray values in 
the range from 0 to 255). Since we divide by g in right-hand side of the Monge-Ampere 
equation (3.7), the function g should be bounded away from zero. To guarantee this 
lower bound we adjust the image and use the modified function 


( 6 . 1 ) 


5 (Z) := 5 (Z) -h max{0, L - min g{Z')} 


with L G N, see also O (5.9)]. Numerical experiments show that the value L = 20 
leads to good results. In order to satisfy the energy conservation condition (3.2) the 
function g needs to be scaled accordingly. 


7. Simulation results. In this section we discuss some numerical simulation 
results obtained by the collocation method for the inverse reflector and refractor 
problems. 
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e 


(a) Luminous intensity as a function on the 
angle 0 between z-axis and direction of ob¬ 
servation. 



(b) Intensity as density 
function / : > R+. 


Figure 7.1: Light emitting characteristics of the radiator of Lambertian type. 



Figure 7.2: Geometrical setting of the examples for the reflector problem. 


7.1. Lambertian radiator and target illumination. For both optical prob¬ 
lems and all of our simulations we use the domain [/ = {x S 5^ : x' S (— 
and a light source with a Lambertian-type emission characteristics. Its emitted lumi¬ 
nous intensity 1(9) is rotationally symmetric, shows a fast decay and is proportional 
to cos(^0), where 9 S [0, JjTt] is the angle between the z-axis and the direction of 
observation. Figure [Td] shows the emission density function / depending on our two- 
dimensional parameter x' and on 9. We choose a light source with this characteristic 
because the maximum possible angular direction for our rectangular domain is about 
^*max = 25° and we therefore have a very low intensity at the edges of fl to make the 
setting more difficult. 

As desired target illumination patterns we chose four images with a variety of 
characteristics, i.e. many different patterns and features, see first row in Figure 7.4 


The first three test images are taken from [9], while the fourth test image is our 
institute’s logo. 


7.2. Geometrical setting and verification. Figure [7i2| shows our geometrical 
setting for the inverse reflector problem where the resulting reflectors are approxi¬ 
mately of the size as in this figure. Here we have S = [4, 12] x [—4,4] x {20}. 

For the refractor problem the dimensions including the size of the optical surfaces 
are chosen very similarly to the case of the reflector problem to have a comparable 
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Figure 7.3: Geometrical setting of the examples for the refractor problem. 


situation, see Figure [773| Here we use a part of the surface of a sphere with radius 0.5 
as initial guess, see also Section 6.6 and S = [—4,4] x [—4,4] x {20}. As refractive 


indices we use ni = | for the lens representing an average glass material and 712 = 1 
for the environment. 

The calculated reflector or lens is verified using the ray tracing software POV- 
Ray [7]. 

7.3. Choice of the parameters. In the nested iteration we successively solve 
the nonlinear systems of equations for the following pairs {N, n) of grid resolutions: 
(16,163), (31,163), (31,55), (61,55), (61,19), (121,19), (121,7), (241,7), (241,3), and 
(481,3), see Section 6.5 for the details. The Newton-type method ends after at most 
200 iterations. 


The regularization parameter in the modified determinant as defined in Sect ion [6f2 
is set to A = 10^ which turns out to be an appropriate choice for all examples. 


7.4. Results. The results of the numerical simulations are depicted in Fig¬ 


ure 


7.4 In the first row the original test images are shown. The first three of them are 


chosen to examine different characteristics within the images, like thin straight lines 


and lettering as in the image “Boat”, see Figure 7.4 (a) Different patterns of high and 
low contrast are present in the image “Goldhill”, see Figure 7.4 (b) in particular at the 


windows and roofs of the houses and the surrounding landscape in the background. 
The image “Mandrill” in Figure 7.4 (c) shows the face of a monkey with a lot of fine 


details like the whiskers. The fourth and most challenging of our test pictures is the 


logo of our institute in Figure 7.4 (d) because it shows the highest possible contrast 


and contains jumps in the gray value from black to white. The iteration counts and 
timings for the numerical experiments are given in Tables [TT] and [7?^ 

First, we notice that for a given original image the output images obtained by 
forward simulation for the reflector and refractor problem look very similar. In com¬ 
parison to the original images the output images are slightly blurred and have a little 
less contrast but visually they only differ locally at very few locations. Major devia¬ 
tions can be observed in the background of the institute’s logo which is not completely 
black after the forward simulation of the mirror and the lens. This is because of the 
minimal gray value needed to avoid the division by zero, see Section [6.7[ 

We see that all of these characteristics of the first three test images are well 
preserved by our method. The computing time for the refractor is approximately 
twice as long as for the reflector but still acceptable with about 4 minutes. 

For the fourth test image we had to adjust the parameters in the nested iteration 
process to handle the sharp edges and work with a finer grid, see Table [7?^ We also 


raised the minimal gray value in Section 6.7 from 20 to 30 obtaining a proportion 
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(a) “Boat” 


(b) “Goldhill” (c) “Mandrill” (d) Institute’s logo 


Figure 7.4: Simulation results for three test images for the reflector and refractor 
problem, first row: desired distribution (original image, image sizes are 
512 X 512 pixel for the first three and 988 x 988 pixel for the last image); 
second row: distribution after forward simulation by ray tracing for the 
reflector problem (result); third row: same as second row but for the 
refractor problem (result). 


between black and white of 1 : 9.5. These parameters lead to results showing also 
a very sharp logo for both the inverse refractor and reflector problems. Note that 
nevertheless in two stages of the nested iteration for the refractor problem the quasi- 
Newton method was stopped because the maximal number of iterations was reached 
without meeting the required tolerances, see Table 7.2 This happened only for two 


intermediate steps of the nested iteration process while we observe convergence in the 
last iteration, which shows us that this does not affect the overall method. In the case 
of the refractor problem the gray line below the letters is irregularly illuminated and 
slightly too bright. Nevertheless, the shape of this line is reproduced very precisely. 


The optically active surface of the lens for the projection of the institute’s logo is 
displayed in Figure [73] Note that the characters used in the logo can be recognized 
on the surface. We observe that they cover about the half of the lens’ surface while 
this is not the case in the original image. Of course this is what we expect because 
we want to redirect a maximal amount of incoming light onto these letters. 
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Table 7.1: Number of iterations of the Newton-type method for each of the ten nested 
iterations and overall computing time in seconds for the standard test im¬ 


ages in Figure 7.4 (a) 
” Mandrill “. 


Boat“, Figure 7.4(b) ”Goldhill“, and Figure 7.4(c) 


Iterations 

{N,n) 


refractor 



reflector 


Boat 

Goldhill 

Mandrill 

Boat 

Goldhill 

Mandrill 

(16,163) 

70 

65 

79 

13 

13 

11 

(31,163) 

13 

13 

15 

13 

11 

11 

(31, 55) 

24 

15 

15 

13 

13 

13 

(61, 55) 

43 

15 

15 

13 

13 

13 

(61, 19) 

35 

44 

37 

13 

13 

13 

(121, 19) 

41 

32 

43 

13 

13 

13 

(121, 7) 

47 

41 

42 

13 

18 

18 

(241, 7) 

46 

38 

40 

13 

15 

13 

(241, 3) 

39 

41 

54 

15 

15 

26 

(481, 3) 

38 

36 

45 

15 

15 

24 

Time / s 

227 

210 

234 

90 

90 

133 


Table 7.2: 


Number of iterations of the Newton-type method for 
iterations and overall computing time in seconds for 


Figure 7.4 (d) 


each of the ten nested 
the institute’s logo in 


Iterations 
{N, n) 

refractor 

reflector 

Institute’s logo 

Institute’s logo 

(21,100) 

33 

54 

(41,100) 

13 

11 

(41,100) 

11 

11 

(81,100) 

13 

11 

(81, 73) 

54 

19 

(161, 73) 

35 

13 

(161, 25) 

200 

90 

(321, 25) 

66 

20 

(321, 9) 

200 

155 

(641, 9) 

59 

22 

Time / s 

1390 

928 



(a) Refractor surface in correct geometri¬ 
cal position (overview) 


(b) High-frequency components of the re¬ 
fractor (fine structure). 


Figure 7.5: Outer refractor surface for projecting our institute’s logo. 























DESIGNING ILLUMINATION LENSES AND MIRRORS 


15 


8 . Summary and outlook. For the efficient and stable solution of the inverse 
reflector and refractor problems we propose a numerical B-spline collocation method 
which is applied to the formulation of the inverse optical problems as partial differen¬ 
tial equations of Monge-Ampere type and appropriate boundary conditions. Several 
challenges for the construction of a stable numerical solution method have been met, 
e.g. we detailed how to enforce ellipticity constraints to ensure uniqueness of the 
solution and how to handle the involved boundary conditions. A nested iteration 
approach simultaneously considerably improves the convergence behavior and speeds 
up the numerical procedure. 

For the inverse refractor problem our algorithm provides a reliable and fast 
method to compute one of the two surfaces of the lens under the assumption of a 
point-shaped light source. Shaping the second surface of the lens, e.g. to minimize 
Fresnel losses, and exploring possible solution strategies for the problem for extended 
real light sources are topics of upcoming research. 

Acknowledgments. The authors are deeply indebted to Professor Dr. Wolfgang 
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of Monge-Ampere type. We thank Elisa Friebel, Silke Glas, and Gudula Kammer for 
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